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Supersonic Reacting Internal Flow Fields 

J. Philip Drummond* 


The national program to develop a trans-atmospheric vehicle has kindled a renewed 
interest in the modeling of supersonic reacting flows, both in the United States and abroad. 
A supersonic combustion ramjet, or scramjet, has been proposed to provide the propulsion 
system for this vehicle. Work has been underway for the past 25 years to develop and 
optimize a scramjet propulsion system for a variety of purposes, but during most of this 
period, the program has not reached the level of intensity that it currently enjoys. With the 
maturing of the scramjet development program, techniques to model the engine flow field 
have also progressed significantly. This has been due both to an advancement in numerical 
methods for computing reacting flow fields as well as an appreciable growth in computer 
speed and storage. The improvements in both algorithms and computer power have also 
lead to a gradually improved understanding of the physics of reacting flow fields which is 
so very important if computation is to have a truly significant impact on scramjet design. 
This chapter will deal principally with the development of computational techniques for 
modeling supersonic reacting flow fields, and the application of these techniques to an 
increasingly difficult set of combustion problems. Since the scramjet problem has been 
largely responsible for motivating this computational work, we will begin with a brief 
history of hypersonic vehicles and their propulsion systems. This will be followed by 
a discussion of some early modeling efforts applied to high speed reacting flows. We will 
then move to the present day and discuss current activities to develop accurate and efficient 
algorithms and improved physical models for modeling supersonic combustion. Following 
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that discussion, some new problems where computer codes based on these algorithms 
and models are being applied will be described. We will then be ready to look beyond 
the current day and draw some conclusions concerning future needs and directions for 
modeling supersonic reacting flows, and hopefully challenge the reader to tackle one of 
these exciting problems that we face in the future. 

Introduction 

Research to develop a supersonic combustion ramjet, or scramjet, propulsion system 
was underway in the late 1950’s. In unrelated efforts, work had also begun to develop 
computational techniques for solving the equations governing the flow through a scramjet 
engine. The marriage of scramjet technology and computational methods for assisting in 
its evolution would remain apart for another decade, however. The principle barrier to the 
union was the lack of a high-speed computer technology for solving the discrete equations 
provided by the numerical methods. Computer resources remain even today as a major 
pacing item in overcoming this barrier. Significant advancement has been made over the 
past thirty years, however, to model the supersonic chemically reacting flow in a scramjet 
combustor. To see how the the two fields finally merged, it is useful to briefly trace the 
evolution of the technology in both areas. 

Following a moderate level of activity in the late 1950*8, there was a significant in- 
crease in the research to develop scramjet engine concepts in the 1960’s. In 1965, the 
NASA Langley Research Center initiated the Hypersonic Research Engine (HRE) Project 
to develop a high speed airbreathing technology base that could then be applied to the 
development of propulsion systems for hypersonic cruise vehicles [l]. The goal of the HRE 
Project was to flight test a regeneratively cooled, hydrogen fueled, pilon mounted scramjet 
on the X-15 research airplane and demonstrate design performance levels. The HRE did 
not reach the flight demonstration stage due to cancellation of the X-15 program, but the 

ground based program did continue and resulted in the development and construction of 
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two variable geometry engine models. Work with these models significantly increased the 
scramjet technology base to be applied in more advanced configurations. 

Following completion of the HRE Project, attention moved to propulsion concepts 
that would provide high performance when installed on a vehicle. The pilon mounted 
HRE would have resulted in excessive levels of external drag, so the pilon was removed 
and work began to highly integrate the engine with the airframe of candidate vehicles. 
In addition, engine weight was reduced by moving from a variable to a fixed geometry 
which reduced the amount of engine structure. Out of this activity, the Langley airframe 
integrated scramjet engine concept was conceived. That program, which has continued to 
the present day, resulted in the successful demonstration of the concept to produce net 
thrust in subscale hardware. A detailed review of this program was given by Northam et 
al. in reference [l], and the reader is referred to their discussion for further details. 

In addition to the NASA scramjet research and development program, other govern- 
ment activities included a Navy sponsored scramjet program at the Applied Physics Lab- 
oratory of the Johns Hopkins University (JHU/APL) (2,3]. This work also increased in 
the 1960’s and was directed towards the development of an air-breathing shipboard missile 
utilizing a scramjet propulsion system. Development of this concept continued until 1977. 
At that time, concern over the storage of highly reactive and toxic fuels to be used by the 
system forced a change to more conventional but safer fuels. This change resulted in the 
development of an integral rocket/duel combustor ramjet concept that utilized a fuel-rich 
gas generator to preburn the fuel for a main supersonic combustor, thus allowing the use 
of hydrocarbon fuels [4], 

The Air Force also sponsored scramjet research and development during the 1960’s [2], 
They continued the support of several programs that were initially funded by the HRE 
program. In 1964, a program was started at the General Applied Science Laboratory to 
continue development of a low-speed fixed geometry scramjet engine, and a duel-mode 
scramjet program was continued with the Marquardt Company at the same time. Soon 
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thereafter in 1965, the Air Force began an effort with the United Aircraft Research Labora- 
tory to continue development of a water cooled variable geometry scramjet design. These 
three efforts ended in 1968, and only the NASA and JHU/APL programs continued into 
the 1970’s. 

It wa s also during the 1970’s that computational techniques were first applied to study 
the supersonic reacting flow found in a scramjet combustor. A detailed review of those 
activities was given by White, et al. in reference [2], and a summary of that discussion 
and additional work is now provided. Some of the earliest work to model supersonic re- 
acting flows was undertaken by Ferri [5] and his colleagues, Morretti [6], Elelman [7], 
and Dash [8,9]. They employed an explicit viscous characteristics method that split the 
governing equations into hyperbolic and parabolic parts followed by a coupled numerical 
solution of each part at each integration step. Modeling of multistep finite rate chemistry 
was also included in their solution strategy. Spalding and his colleagues then took Ferri’s 
splitting-based approach and improved its efficiency by developing a fully implicit solu- 
tion procedure for solving the governing equations [10]. Spalding then developed several 
implicit parabolized Navier-Stokes programs for modeling scramjet combustor flow fields. 
These codes included the CHARNAL two-dimensional axisymmetric program [11] and 
the SHIP three-dimensional program [12]. Both programs utilized the well known SIM- 
PLE solution procedure for spatially marching the governing equations in the parabolized 
direction while employing a tri-diagonal matrix solution procedure to perform repetitive 
sweeps for solution of the equations in the cross-plane(s) [13]. The programs were written 
to assume that a state of chemical equilibrium alway existed, but they were later modified 
by Evans [14] to include the effects of finite rate chemical reactions. The modified pro- 
grams are still being used today for studies of mixing and reaction in candidate combustor 
configurations. 

The work of Ferri and Spalding was then adapted by Dash to develop the SCORCH 
program that used a hybred explicit/implicit procedure for modeling supersonic reacting 
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flows. The method again split the governing equations into hyperbolic and parabolic parts. 
The hyperbolic part was solved using a viscous characteristics approach that employed 
an upwind finite difference procedure. The parabolic part was solved using an implicit 
finite difference procedure [15]. Work on this program and its application to supersonic 
combustion problems has continued to the present day. 

While Ferri, his colleagues, and Spalding were developing analysis techniques for di- 
rect application to the supersonic reacting flow problem in a scramjet, other algorithm 
development work was underway, directed primarily at solving high speed external flow 
problems. These techniques ultimately found their way, however, into the internal reacting 
flow arena. The first of these algorithms was the MacCormack explicit, unsplit predictor- 
corrector method that was initially developed to model the hypervelocity impact cratering 
problem [16]. The MacCormack method was a variation of the Lax- Wendroff second order 
accurate scheme. The method was robust and easily applied to complex geometries. Be- 
cause of these qualities, the algorithm was readily adopted and used to study a wide class 
of external flow problems. In fact, due to its general applicability, MacCormack ’s unsplit 
algorithm is still used today as one option in several codes that are applied extensively to 
the modeling of scramjet flow fields. Implicit algorithms were also developed for external 
flow problems in the 1970’s. Their development was motivated by the need to resolve the 
high gradients present in wall boundary layers. The resolution of boundary layers required 
fine computational grids, resulting in a severe stability constraint on the marching step 
size of an explicit method. Where only a steady state solution was required, i.e., time 
accuracy was not necessary, implicit methods could achieve a significantly more rapid rate 
of convergence. Early work to develop implicit solution techniques for the Navier-Stokes 
equations was carried out by Briley and McDonald [17] and Beam and Warming [18], 
Both approaches used a spatial factoring procedure that reduced the multidimensional 
problem to one of sequentially solving a set of one dimensional spatial implicit operators. 
Using this computationally efficient procedure, convergence rates one to two orders of mag- 
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nitude faster than the explicit method were achieved for steady-state problems on highly 
stretched grids. 

While the application of implicit methods was generally limited to scramjet inlet flow 
fields through the late 1970’s and early 1980’s, explicit methods were applied extensively in 
studies of combustor flow fields. The MacCormack method was employed by Drummond to 
model internal scramjet combustor flow fields. In 1977, he developed the two-dimensional 
TWODLE combustion program based on that method. The code used an equilibrium 
chemistry scheme to model H 2 -air reaction and several algebraic eddy viscosity methods 
to model the turbulence field. The program was applied to several scramjet combustor 
component problems. Particular emphasis was given to the scramjet fuel injector problem 
in an attempt to better understand the complex flow field in this region of the engine 
[19,20]. Development on the program continued into the early 1980’s when the program 
was used to carry out the first simulation of a scramjet flow field using a two-dimensional 
model engine module [22]. Detailed studies to optimize candidate scramjet fuel injector 
configurations were also completed during this period [21,23]. 

An explicit solution procedure was also employed by Schetz during the early 1980’s to 
model the APL duel combustion ramjet described earlier [26], He employed a modular 
approach to carry out his analysis. The mixing and burning of the center jet from the 
fuel-rich gas generator was calculated with a jet mixing code [24,25] that was modified 
to include a turbulent kinetic energy turbulence model, a chemistry model, and other 
improvements. Because of the high static pressures and temperatures that were present 
in the device, a local diffusion-controlled, equilibrium chemistry model was used to model 
reaction in the combustor. Schetz’s procedure for modeling combustor flows was ultimately 
combined with an inlet analysis procedure to compute performance estimates for the duel 
combustion ramjet [27]. 

While numerical methods for modeling scramjet flow fields were developing through 

the 1960’s, 1970’s and early 1980’s, there was a parallel growth in computer hardware 
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upon which these methods could be applied. Many of the early calculations were carried 
out on IBM 7090 and CDC 6600 class machines. Hardware improvements, that allowed 
the consideration of more realistic problems came in the late 1960’s with the arrival of 
the CDC 7600 computer. The most significant hardware improvement came in the mid 
and late 1970’s, however, when vector processing supercomputers became available to the 
computational community. These machines included the CDC Star-100 and the Cray 1. 
They were followed in the early 1980’s by the Cyber 205 and the Cray X-MP which gave 
performance capabilities several orders of magnitude greater than the scalar machines 
available to the researcher prior to their arrival [2]. To this time, the state of computer 
resources had resulted in a major barrier to advancing the state of the art in modeling 
supersonic reacting flows. With the new machines, however, the researcher was now in 
a position to begin dealing with the detailed physics contained in these complex flows. 
The burden now returned partially to the state of numerical algorithms used to model 
supersonic combustion. This state has continued to the present day. 

We are now in a situation where both numerical algorithms and computer technology 
are pacing our ability to formulate an improved understanding of supersonic reacting flows. 
We will concentrate for the remainder of this chapter on the numerical challenge, i.e. what 
is needed to advance the computational state of the art that will result in an improved 
understanding of supersonic reacting flows. We will then explore how we can use this 
improved understanding to solve practical problems associated with the modeling and 
design of a scramjet combustor. There is a critical need today for creditable methods for 
modeling flows typical of those found in a supersonic combustor. The National Aero-Space 
Plane, mentioned earlier in this chapter, will operate at Mach numbers as high as 25. To 
develop a successful design for the propulsion system of this vehicle, extensive use will be 
made of ground-based facilities to create flow fields consistent with those that the engine 
will injest over its operating envelope. Unfortunately, however, ground based facilities 
are only able to create continuous flow conditions up to a flight Mach number of about 
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8. Beyond Mach 8, the only options available to the experimentalist are pulse facitities 
that create flight conditions for only a short period of time, providing a data collection 
window of only a few milliseconds. Numerical methods provide an alternative to the Mach 
8 barrier, but only if they are properly applied to the problem. To examine the challenge 
that this poses for those who are applying computational methods, we will proceed along 
the following path. We will first review the equations that govern the supersonic reacting 
flow problem and the modeling that these equations require. We will next explore a number 
of promising numerical methods, both old and new, for accurately solving these governing 
equations. Several solutions for reacting flow problems using some of these methods will 
then be presented to assess the capabilities of the techniques and the computers which 
provided their results. We will then be in a position to evaluate where we are with these 
methods and where we need to go. That will be the subject of the conclusion to this 
chapter, or at least this author’s perspective of it! 

Theory 


Governing Equations 

The Navier-Stokes, energy, and species continuity equations governing multiple species 
undergoing chemical reaction have been derived by Williams [28] . The terms used in 
these and subsequent equations are defined in the appendix at the end of this chapter. 
The governing equations are given by 

Continuity 
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Energy 
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Species Continuity 
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Radiation heat transfer is not included in equation (6). Also, 
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The diffusion velocities are found by solving 
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Note that if there are ns chemical species, then i = l,2,...,(ns-l) and (ns-1) equations must 
be solved for the species / f . The final species mass fraction /„, can then be found by 
conservation of mass since YJiLi /,* = !. 
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Thermodynamics Model 


To calculate the required thermodynamic quantities, the specific heat for each species 
is first defined by a fourth-order polynomial in temperature 

= Ai + BiT + CiT 2 + AT 3 + EiT 4 (11) 

The coefficients are found by a curve fit of the data tabulated in reference [29]. Knowing 
the specific heat of each species, the enthalpy of each species is then found from equation 
(8) and the total internal energy is computed from equation (7). 

To determine the equilibrium constant (required in the next section) for each chemical 
reaction being considered, the Gibbs energy of each species must first be found. For a 
constant pressure process, ^ from equation (11) is first integrated over temperature to 
define the entropy of the species, and then the resulting expression is integrated again 
over temperature to obtain a fifth-order polynomial in temperature for the Gibbs energy 
of each species. 


a 
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MT-T In T) + § T 1 + ~T 3 + 2L t * + %-T* + F, - G<T 


12 


20 


( 12 ) 


The coefficients F{ and G{ are again defined in reference [29], The Gibbs energy of 
reaction is then calculated as the difference between the Gibbs energy of product and 
reactant species. 


ng n» 

AGn y = J2 Si - Y. liiSi j = 1, 2, nr (13) 

«=1 *=1 

The equilibrium constant for each reaction can then be found from [30] 

i —AC 1 

in) 

where An is the change in the number of moles when going from reactants to products. 
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Chemistry Models 


The rate of reaction of chemical reactions is often defined by using the Arrhenius Law. 
A modified form of of the Arrhenius law is usually employed when modeling supersonic 
combustion. It is given by 


K fi = AjT N > exp( 


z El 

T 


(15) 


The values of the preexponential constant A, power constant N, and the activation energy 
E have been determined for a number of reaction schemes. Unfortunately, there is a great 
deal of uncertainty for many chemical reactions. One of the best understood mechanisms, 
however, is the hydrogen-air reaction system. This is not the reason that hydrogen fuel 
was chosen for several scramjet concepts, but it has proven convenient for its combustor 
analysists! Values for A, N, and E for a typical hydrogen-air mechanism are given in Table 
1. Knowing the forward rate, the reverse rate is then given by 


K\,. 




K 


«<?/ 


(16) 


Once the forward and reverse reaction rates have been determined, the production rates 
of the species are found from the law of mass action. For the general chemical reaction 


ns Kt. ns 

<=1 i=l 


j ~ 1,2, ..., nr 


(17) 


the law of mass action states that the rate of change of concentration of species i by reaction 
j is given by [28] 


(C'.-)i = hi, - %•.) [K fi n <?>■ - K hj n cl * ] i = 1,2,..., ns (18) 

*=1 1=1 

The net rate of change in concentration of species i by reaction j is then found by summing 
the contributions from each reaction 


c, = £(Q),- 

i-i 


( 19 ) 
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Finally, the production of species i can be found by multiplying its rate of change of 
concentration by its molecular weight. 

w { - CiMi ( 20 ) 

The source terms in equation (4) are now determined as a function of the dependent 
variables. 


Molecular Diffusion Models 


The coefficients governing the molecular diffusion of momentum, energy, and mass are 
determined from models based on kinetic theory. The set of models that is often used is 
now described. Individual species viscosities are computed form Sutherland’s law 


Vo K T 0 } T+S 


( 21 ) 


where /x 0 and T 0 are reference values and S is Sutherland’s constant. These constants are 
tabulated for many species in references [31,32]. Once the viscosity of each species has 
been determined, the mixture viscosity is found from Wilke’s law [33] 


Vm 


1=1 


Vi 


+ i Xjfri 


( 22 ) 


where 


<t>u = i 


1 + (£“> M (#) M T 


(23) 


Species thermal conductivities are also computed from Sutherlands’s law 


- = ( T ) i* T ° + S ' 
K k tJ T + S' 


(24) 


with different values of the reference values k 0 and T' c and the Sutherland’s constant S' . 
These values are also tabulated for a number of species in references [31,32]. The mixture 
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thermal conductivity is computed using conductivity values for the individual species and 
Wassilewa’s formula [34] 


ki 


km — Yl 7 . 1 y*»/ x i • 

where <^ iy = 1.065 </>ij and faj is taken from equation (18). 


(25) 


For dilute gases, Chapman and Cowling used kinetic theory to derive the following 
expression for the binary diffusion coefficient j Dy between species i and j [31], 

o.ooisssr 1 -^^^^) 0 - 5 


Ay = 


pcr/ y n D 


(26) 


Here, the diffusion collision integral Hu is approximated by 


n D 


/fi-0.145 


+ (r + 0.5) 


—2 


(27) 


where T — Values of the effective temperature T e and the effective collision diameter 
*•/ 

o are taken to be averages of the separate molecular properties of each species, giving 


Oij - 0.5 (cq 4- o j) 


(28) 


and 


r... = (T'.t ( .) 


0.6 


(29) 


For most molecules, the thermal diffusion coefficient is generally small when compared 
with the binary diffusion coefficient, and therefore, the thermal diffusion coefficient can 
be neglected. This is a fortunate fact, since values of the thermal diffusion coefficient 
are generally not known for most species. For low molecular weight molecules such as 
hydrogen, though, the thermal diffusion coefficient can be important. A set of relationships 
for the thermal diffusion coefficient of species having a molecular weight less than 5 has 
been developed by Kee et al. [35], The reader is referred to reference [35] for further 
information and some numerical details for computing thermal diffusion coefficients of light 
molecules. 
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Once the binary and thermal diffusion coefficients for all species combinations are 
known, the diffusion velocities of each species can be computed from equation (10). The 
diffusion velocity is the velocity induced upon each species by all diffusion processes that are 
present in the flow. The solution of equation (10) requires solving a simultaneous equation 
system, with the number of equations equivalent to the number of species present for each 
component of the diffusion velocity. It should be noted that for i species, however, the 
system of 1 equations defined by (10) is not linearly independent. One of the equations must 
be replaced by the constraint pfiVi = 0 to make the system linearly independent. The 
resulting simultaneous system of equations must then be solved for the diffusion velocities. 

The process of solving for the diffusion velocities can be computationally quite expen- 
sive, A coupled system of equations must be solved for each of the three components of 
the diffusion velocity at each computational grid point. This process can require as much 
time as solving the Navier-Stokes equations for the three components of the convection 
velocities. Alternately, for hydrogen-air chemistry where large amounts of nitrogen are 
present, it is sometimes assumed that each species is present as a ” trace” in a mixture 
with N 2 [26]. Then each species is assumed to diffuse only into N 2 with that process 
defined by its binary diffusion coefficient with N 2 . Finally, for engineering calculations, it 
is often further assumed that the diffusivities of each chemical species present in the flow 
are the same. Then the diffusion of each species into the remaining species varies only with 
its respective concentration gradient. The diffusion velocities then decouple, and equation 
(10) reduces to 


Vi,i = 


Ddfi 
fi dxj 


(30) 


where V u is the diffusion velocity vector of the i th species in the j th coordinate direction 
(j = [x,y,z]) and D is the binary diffusion coefficient. If the binary diffusion with N 2 is not 
used, the value of D is determined by chosing an appropriate value of the Schmidt number 
Sc since D = The mixture viscosity n is determined as before from Wilke’s law. When 

the binary diffusion assumption is invoked, it is often further assumed that the mixture 
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thermal conductivity can be defined by k = after an appropriate value of the Prandtl 
number Pr has been chosen. 

Turbulent Diffusion Models 

While the techniques for defining the molecular diffusion of momentum, heat, and mass 
are reasonably well established in a supersonic reacting flow, the same statement cannot 
be made for our ability to describe the turbulent diffusion of these quantities. Work to 
develop methods for modeling turbulent supersonic combustion is now in its early stages. 
Conventional approaches have included the use of algebraic eddy viscosity models or dif- 
ferential transport models. Several eddy viscosity models have been used, in particular the 
Cebeci-Smith model [36] and the Baldwin-Lomax model [37]. The differential transport 
models include the k/c turbulent kinetic energy model and its variants [38], a modified 
k/e model that included a supersonic flow compressibility correction [39,41], and a multi- 
ple dissipation length scale (k/multiple c) model with a compressibility correction [39,40] 
that addressed the existence of multiple dissipation length scales that exist in the energy 
cascade of a turbulent flow. In addition to these differential transport models, the alge- 
braic Reynolds stress models of Rodi [42] and Sindir [43] have also been considered for 
use in modeling turbulent supersonic reacting flows. A review of all of these models has 
been given by Sindir [43]. In that review, he also critically compared the models against 
several nonreacting flow experiments prior to using the models for studying flows with 
reaction. He concluded that forms of the algebraic Reynolds stress model that he consid- 
ered produced the best agreement with nonreacting data. He also found that the multiple 
dissipation length scale model did not offer any advantage over the basic k/e model. 

All of the turbulence models described above have a major disadvantage when applied 
to reacting flow fields. They fail to account for the important coupling between the fluid 
mechanics and the chemistry. Turbulent fluctuations in the fluid mechanic variables have 
a direct effect upon the species production rates. The coupling between these two fields 
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occurs through the Arrhenius rate expression, equation (15), and the law of mass action, 
equation (18). The Reynolds averaging process applied to the governing equations elimi- 
nates the direct effect of temperature and species fluctuations on species production rates. 
For example, a positive temperature fluctuation would cause a decrease in the size of the 
exponential argument of the Arrhenius rate expression, with a corresponding increase in 
the forward kinetic rate of a particular reaction. This would in turn produce an increase in 
the time rate of change of the products of that reaction. More importantly, if the reaction 
were at a critical stage, where perhaps a small increase in temperature would cause a reac- 
tion to enter an ignition stage, the entire species distribution of the flow field downstream 
could be changed. 

Two promising ways for accounting for the effects of fluid and species fluctuations on 
chemical reaction would be through probability density functions or direct numerical sim- 
ulation. The application of the probability density function approach to a reacting flow 
has been covered by Stephen Pope in a companion chapter in this book and so that subject 
will not be further addressed in this chapter. Direct numerical simulation offers another 
attractive approach for modeling a turbulent reacting flow. The method has been used 
for a several years to accurately model lower speed reacting flows [44,45,46]. With this 
approach, the Navier-Stokes and species continuity equations are resolved down to the 
smallest scale features of the flow field. The size of those scales goes inversely with the 
Reynolds number of the flow field. Clearly then, for the high Reynolds numbers that occur 
in typical supersonic reacting flows, the smallest scales can become quite small, nessitating 
a very fine computational grid to resolve the scales. Also, when high speed flow undergoes 
chemical reaction, additional scales are introduced by the combustion process. Herein lies 
the principal difficulty of applying direct simulation to a high speed flow. The difficulty is 
not so much one of numerical algorithms as it is of computer power. Highly accurate numer- 
ical algorithms are required, but appropriate high-order finite-difference/volume methods 
or spectral methods have been developed that satisfy that requirement. The large number 
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of computational grid points required to resolve the smallest scales in the flow requires 
large computer storage, and therefore, meaningful calculations can be carried out only- 
on large memory machines. Currently, direct numerical simulations have been made for 
nonreacting flows with Reynolds numbers up to about 10000 on a Cray 2 computer [47]. 
Work is proceeding to directly simulate a chemically reacting flow of a similar Reynolds 
number, and those activities will be discussed with other applications later in this chapter. 

As an alternative to direct numerical simulation with its intensive memory require- 
ments, it is possible to model rather than compute the smallest scales. In this approach, 
termed large eddy simulation, the larger scales above a chosen wavelength are still com- 
puted. The smaller scales below the cutoff wavelength are modeled, however, using a 
subgrid scale model. Large eddy simulation is an attractive alternative because only the 
larger scale effects are computed, lessening the computer memory requirements for higher 
Reynolds number flows. Subgrid scale models must be constructed, though, that give an 
accurate rendering of the physics of small scale phenomena. This is a difficult task. Work 
is underway to develop subgrid scale models for nonreacting flow, for example the early 
work of Schumann [48] and later work described by Speziale et al. [49]. Large eddy 
simulation is an attractive technique for modeling high speed reacting flows. Little has 
been done so far with this technique, but it warrants serious attention in the future. 

Discretization of the Governing Equations 

Once the governing equations and required modeling are in hand, the numerical method 
of choice can be applied to discretize the governing equations in space and time. The 
numericist has three basic options for discretizing the equations in time. He may express 
the equations explicitly, implicitly, or in a partially implicit manner. The merits of the first 
two approaches were discussed in the introduction. The latter approach is attractive when 
the time scales for chemical reaction are quite small as compared to the prevailing fluid 

dynamic time scales. In this case, the governing equations become stiff, and a significant 
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advantage can be gained in convergence of the equations to steady state by casting only 
the source term in the equations implicitly [50,51]. Before discretization, it is convenient 
to express the governing equations (1) in vector form. In that form they become, 

dO dP dP dG 


-f- • -4- ' -4- — — — JJ 

dt dx dy dz 


(31) 


where U is the vector of dependent variables, E y F, and G are flux vectors continuing 
convective and diffusive terms, and B[ is the source term containing body forces and 
the chemistry production terms. The temporally discrete form of equation (31), written 
explicitly, is then given by 


ff» + i = + 


dE n dF n dG n 
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where n is the old time level and n+1 is the new time level. Written implicitly, equation (31) 
becomes 


— + ^ + 1 (33) 

A partial implicit statement of equation (31) is obtained when only the source term is 
written implicitly, i.e. JT n+1 , and the remaining terms are written explicitly. 

Once the temporal discretization of equation (31) has been chosen, the spatial deriva- 
tives must also be discretized. There are many choices available. In the next section of this 
chapter, we will examine a number of those choices, using earlier techniques as well as some 
newer ones. We will then be in a position to review some applications of these methods to 
practical supersonic combustion problems that we are able to numerically simulate today. 


Numerical Algorithms 


A number of numerical algorithms have been used over the past 20 years to solve the 
equations that govern a supersonic reacting flow field. The earliest of those approaches 
were reviewed briefly in the introduction to this chapter. We will now discuss, in somewhat 
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more detail, several of those algorithms that are still in use today. New algorithms that 
have appeared fairly recently will then be considered, and their merits and the advantages 
that they offer over the older approaches will be discussed. Accurate methods used in 
other fields and recently borrowed to model combustion problems will also be reviewed. 
With the review behind us, we will then move on in the next section to the application of 
these algorithms to several practical combustion problems. We will then be in a position 
to assess where we are in our ability to model these practical combustion problems and 
what is needed to further extend our capabilities. 

Conventional Approaches 

The algorithms developed by Spalding, Dash, MacCormack, and their colleagues today 
continue to be popular tools for modeling supersonic reacting flows typical of those found 
in scramjet combustors. Each of these algorithms was described in the introduction, and 
references were given to provide more details. The Spalding three-dimensional parabolized 
Navier-Stokes code, SHIP [12] as modified by Evans [14], is still being used to carry out 
engineering design studies of scramjet configurations as well as basic high speed fuel-air 
mixing studies. The two-dimensional parabolized Navier-Stokes code, SCORCH, of Dash 
[15] has recently seen considerable use to perform analyses of the National Aero-Space 
Plane (NASP) propulsion system. In addition, the SCORCH code has also been used to 
carry out several fundamental studies of experiments being used to design that propulsion 
system. The MacCormack algorithm was employed by Drummond in the TWODLE code 
[19,22] to solve the two-dimensional Navier-Stokes equations describing a scramjet flow 
field. We will examine some interesting computations from these conventional algorithms 
in the applications section of this chapter. 

A number of other extensions of the MacCormack algorithm were made following the 
work that was just described. Drummond extended the TWODLE code [52] to include 

detailed models for finite rate chemistry and kinetic theory based models for the molec- 
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ular diffusion of momentum, heat, and species. He also added an option for treating the 
chemical source term implicitly, as suggested by Bussing and Murman [50j, to allow for a 
more efficient treatment of stiff kinetic source terms. In that form, the new code (SPARK) 
then solved in two-dimensions without simplification the complete governing equation set 
given in equations (1) through (29). Options were also provided, however, to simplify the 
diffusion modeling to use the approach indicated by equation (30). Following development, 
the SPARK code was also applied to NASP configurations, but perhaps more importantly, 
it was also applied to a number of basic high speed reacting flow problems to seek an im- 
proved understanding of important physical processes that occur in these flows and that 
ultimately effect the performance levels that can be achieved by the propulsion system. 

As scramjet technology evolved, a critical need developed for a three-dimensional anal- 
ysis tool for modeling high speed combustor flow fields. Uenishi and Rogers [53] extended 
the three-dimensional nonreacting Navier-Stokes inlet code (NASCRIN) developed by Ku- 
mar [54,55] to include multiple species, but initially they did not include chemical reaction. 
The program again used the unsplit MacCormack method [16] to integrate the governing 
equations (1) through (10). Thermodynamic properties were also defined using equa- 
tion (11) and the procedure described in the discussion following that equation. Molecular 
diffusion of momentum, energy, and species was modeled with Sutherland’s and Wilke’s 
law, the Reynolds analogy, and Fick’s law, respectively. Turbulent diffusion was modeled 
with the Baldwin-Lomax turbulence model [37] along with chosen values of the turbulent 
Prandtl and Schimdt numbers. A clever storage scheme was also employed with Mac- 
Cormack’s predictor-corrector scheme that saved predictor values locally only until the 
corrected values could be computed [55], The scheme allowed computations to be carried 
out, with three chemical species, on computational grids with up to 500,000 points on a 
Cyber-205 having 32 million words of storage. When the code was completed, a number of 
studies were made to model supersonic fuel-air mixing in scramjet like flows [53], Of par- 
ticular interest were the calculations of the downstream mixing of a transverse hydrogen 
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fuel jet injected across a supersonic air flow. These calculations by Uenishi were the first 
simulations of the three-dimensional near-field fuel injector problem in a scramjet engine. 
At that time, a clear understanding of the flow field that existed near the fuel injectors 
was critical to achieving a successful engine design. Work to optimize fuel injector design 
based on these and other simulations is continuing today. 

Uenishi et al. then extended their program to include finite rate chemical reaction 
[56]. Motivated by the need to model hydrogen-air combustion taking place in a scramjet, 
they chose a two-step hydrogen-air reaction model developed by Rogers and Chinitz [57], 
The model considered five species (H 2 , 0 2 , OH, II 2 O y and N 2 [inert]) participating in the 
following chemical reactions. 


H 2 + 0 2 ** 2 OH (34) 

H 2 + 2 OH *=* 2 H 2 0 

The approach defined by equations (15) through (20) was then used with the chemistry 
model to determine values of the chemistry source terms. The source terms were sometimes 
found to be numerically quite stiff, and so the numerical method was also modified to 
include implicit source terms. Otherwise, the numerical method and the physical modeling 
remained unchanged from the previous code. Uenishi again applied his extended program 
to the transverse fuel jet problem, but in this case with chemical reaction. Encouraged by 
those results, he then went on to model an actual combustor configuration. Some results 
from those calculations will be included in the next section of this chapter. 

Because of the need for a more basic modeling capability in three dimensions, Car- 
penter then extended the SPARK combustion code to three dimensions [59]. The three- 
dimensional code retained all of the basic modeling of the two-dimensional code, i.e. it 
solved the system defined by equations (1) through (30). In addition, Carpenter added 
a generalized equilibrium chemistry model and a generalized finite rate chemistry model 
that allowed for consideration of any fuel-air system with any number of reaction paths 
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[60]. Typically, either a seven species, eight reaction model or a nine species, eighteen reac- 
tion model was used to represent hydrogen-air chemistry. The eighteen reaction model is 
given in Table 1. Turbulence was modeled with either a Cebeci-Smith or Baldwin-Lomax 
eddy viscosity model, or a two-equation k/e turbulent kinetic energy model. Calculations 
were then carried out to validate the code by modeling several of the problems considered 
by Uenishi. Following successful agreement with Uenishi’s results and experimental data, 
the code was also used to study a model scramjet combustor. Those results will also be 
presented in the applications section. 


Three-dimensional parabolized Navier-Stokes programs were also developed to model 
supersonic combustor flow fields. These programs often provided a more efficient solution 
procedure if the flow field contained no subsonic regions. The flow field in the neighbor- 
hood of the fuel injectors in a scramjet combustor contains subsonic separated regions 
necessitating a solution of the full (spatially elliptic) Navier-Stokes equations. Somewhat 
downstream of this region, however, the flow takes on a principal supersonic flow direction, 
allowing solution of the parabolized equations and the application of parabolized codes. 
In response to this need, Chitsomboon developed a three-dimensional parabolized Navier- 
Stokes (PNS) program [61] by extending a two dimensional PNS program that he had 
developed earlier [62,63]. He solved the conventional parabolized Navier-Stokes equations 


together with a set of species continuity equations given vectorally by 

d& dP dG 

+ W + ^ = # (35) 

where E, F, G, and H have the same definitions as given in equation (31). The dependent 
variable vector q — [p,pu,pu,ptu,T,/?i, ...] was chosen nonconvention ally, with the energy 
equation written in terms of temperature rather than total enthalpy or total internal 
energy. The equations (35) were then discretized using the Vigneron method [64], and the 
nonlinear vectors F, G, and H were linearized with respect to the dependent variable vector 


q. In order to insure that the numerical scheme was then stable, the flux vector E was 
linearized using the approach of SchifT and Steger [65]. Thermodynamics and chemistry 
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were handled in a manner identical to the approach used in the Uenishi code [58]. Following 
initial completion of the code, it was compared against the Uenishi program and fairly good 
agreement was achieved. The code is undergoing further development today. 

During this same period, Gielda developed a three-dimensional explicit PNS program 
[66] using the MacCormack explicit algorithm [16]. The code was fully vectorized to run 
efficiently on vector supercomputers such as the Cray 2 and the Cyber 205. Gielda found 
that his explicit scheme was quite competitive with implicit algorithms for problems at 
high Mach number or with surface discontinuities. He was able to resolve the problem of 
decoding the axial flux vector (E in equation ( 35)), that had earlier limited the application 
of explicit PNS codes, by also employing the Vigneron procedure of splitting the axial 
pressure gradient. Following completion of this nonreacting program, Gielda extended 
his code (then named the SSCPNS code) by adding the parabolized species continuity 
equations to the governing equation system [67]. He also incorporated both an equilibrium 
and a global one-step if 2 -air finite rate scheme into the program. The extended program 
was then validated against several experimental cases, and generally excellent agreement 
was obtained between data and computation. With this confidence in hand, the code 
was applied to a three-dimensional generic inlet-combustor scramjet configuration that 
included gaseous hydrogen fuel injection and reaction. Some of these interesting results 
will be presented in the following section. 

Kamath then employed the Gielda algorithm to develop a parabolized version of the 
three-dimensional SPARK combustion code [59]. He generalized the coordinate trans- 
formation to allow the streamwise coordinate to be orientated in the most supersonic 
direction. He also utilized the generalized equilibrium and finite rate chemistry schemes 
developed by Carpenter [60] such that any multistep reaction scheme could be considered 
with the algorithm. The extended code was next validated with several test cases, one of 
which was also utilized by Gielda, and it gave generally good agreement with data. Work 
on the program is now continuing. 
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Alternate Approaches 


Several numerical algorithms have been developed to solve the equations governing 
high speed reacting flow fields, but these algorithms have not been applied to the scramjet 
combustor problem. Most of these methods have been developed to model supersonic 
or hypersonic flow with interacting air chemistry moving externally about configurations. 
The remaining algorithms have been developed to study basic phenomena associated with 
high speed reacting flows, but have not yet been applied to scramjet problems. Each of 
the approaches falls into the general class of monotone methods, that is methods that 
employ flux-correcting or flux-limiting procedures to preserve high numerical resolution 
without the numerical oscillations associated with higher accuracy. Included in this class 
of algorithms are Flux Corrected Transport methods, TVD (total variation diminishing) 
methods, and TVD like methods that exhibit TVD behavior. These algorithms would offer 
the modeler advantages over conventional methods when studying scramjet problems, and 
they should be seriously considered for future work. It is for that reason that we discuss 
them here. 

The first monotone method applied to chemically reacting flows was the Flux Corrected 
Transport (FCT) algorithm developed by Boris [68,69,70). In this method, a small amount 
of artificial diffusion is added to the governing equations in smooth regions of the flow to 
stabilize the solution. In regions where high gradients exist, larger amounts of diffusion are 
added to maintain monotonicity. The diffusion is added in such a manner, however, that 
the overall dissipation is held below that of the conventional algorithms. One proceeds as 
follows. Starting with equation (32), the flux terms are discretized in space, and then a 
diffusive term is added to insure positivity as the equations are integrated in time. That 
integration is then performed to determine a first value of the dependent variable vector at 
the next time step. An “anti-diffusive” correction using the first values of the dependent 
variables is then applied to reduce the numerical diffusion added in the first step. Care 
must be taken when applying this step, however, because the anti-diffusive correction can 
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degrade the monotonicity of the method. Therefore, the anti-diffusive terms are limited 
by a “flux correction” procedure such that no new maxima or minima are introduced into 
the solution. The initial fluxes in equation (32) are then replaced by the corrected fluxes 
and the equation is again advanced over the same time step to arrive at the final value of 
the dependent variables at the new time. A more detailed but very readable description 
of the FCT method is given in reference [70]. 

Following development of the FCT method, a more general approach was suggested by 
Zalesak [71]. His approach allowed the method to be readily incorporated into existing 
algorithms that did not provide monotone behavior. In addition, the method could be 
more easily generalized to two and three spatial dimensions. Zalesak viewed FCT as a 
hybridization of a low order and a high order method. The anti-diffusive flux was then 
found as the difference between the fluxes determined by the high order and low order 
methods. Once found, the anti-diffusive flux was then limited as before by flux correction, 
and the solution was advanced using the corrected fluxes. A more detailed discussion of 
Zalesak ’s approach is again given in reference [70]. 

Most of the new work to model high speed reacting flows has taken place over the last 
five years. It was motivated primarily by the need to model hypersonic flows about vehicles, 
including hypersonic cruise aircraft such as the NASP, and reentry vehicles. Therefore, the 
methods were developed to model high speed strongly shocked flows undergoing air chem- 
istry. To compute flows of this type, MacCormack and Candler developed an implicit flux 
split scheme, as an extension to MacCormack’s explicit predictor-corrector finite difference 
method [16], to solve the Navier-Stokes equations. MacCormack initially developed the 
implicit algorithm to consider only nonreacting flows [72]. A finite volume approach was 
used to discretize the flux terms. In addition, Steger- Warming [73] flux vector splitting 
was introduced to more properly account for the propagation of information through the 
flow field. Finally, to relax the constraint on the time step imposed by the Courant con- 
dition in his explicit method, the flux terms were written implicitly and then linearized. 
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This procedure resulted in a coupled set of algebraic equations to be solved at each time 
step once the spatial opperators had been applied to the flux terms. The equation system 
was then solved iteratively using either a line Gauss-Seidel procedure or Newton iteration. 

Following development of the basic algorithm, Candler and MacCormack extended the 
method to consider high speed air flows that were ionized and in thermodynamic and 
chemical nonequilibrium [74,75]. To model such flows, equations describing species con- 
tinuity, vibrational energy of each diatomic molecule, and electron energy were appended 
to the Navier-Stokes equations. The fully coupled system of equations was again solved 
using Gauss-Seidel line relaxation along with an implicit, flux split treatment of the flux 
terms. Air chemistry was modeled with a seven species, six reaction model that included 
^ 2 ,G 2 , NO,NO + , N, O, and e~. For these species, four vibrational temperatures and 
an electron temperature were computed. The chemical source terms associated with the 
seven species and the thermal source terms were also computed implicitly due to the small 
kinetic time scales, relative to fluid scales, that were introduced by the chemistry. When 
the authors completed the extension of their algorithm to include chemical reaction, they 
applied it to several high speed external flow problems in both two- and three-dimensions 
that are described in references [74] and [75]. Even though these examples considered 
only external flows with air chemistry, it appeared that the algorithm could readily be 
modified to consider internal flows with combustion chemistry, and, therefore, serve as a 
means for modeling scramjet combustor flow fields. 

Flux splitting methods were also employed by Grossman, Walters, and Cinnella to 
model high speed chemically reacting flow problems. Grossman and Walters initially de- 
veloped their algorithm to solve the Euler equations for nonreacting flows, but included real 
gas effects [76]. Three forms of flux splitting were considered, including Steger- Warming 
flux vector splitting [73], van Leer flux vector splitting [77], and Roe flux difference split- 
ting [78]. Each of these splitting methods was originally derived to be applied to ideal 
gas flows. They were rederived in reference [76] to allow their application in problems 
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with real gas effects. The new derivations are based on the concept of an equivalent 7 to 
replace the ratio of specific heats for an ideal gas. Each approach is derived with sufficient 
detail in the paper to allow the reader to incorporate the modifications into existing flux 
split codes written to originally model ideal gas flows. Once the flux split equations had 
been developed, they were solved using a two-step predictor-corrector method that was 
second-order accurate in space and time. Spatial differences were formed using the MUSCL 
differencing procedure and flux limiting that was described in reference [79], Following 
the successful application of the algorithm to a one-dimensional shock tube problem, the 
real gas splitting was incorporated into a two-dimensional implicit finite volume code that 
originally utilized van Leer splitting and Gauss-Seidel line relaxation to solve the equations 
governing ideal gas flows [80]. The modified code was then applied to a two-dimensional 
inlet flow field exhibiting equilibrium real gas behavior. 

Grossman and Cinnella then extended the algorithm to include vibrational and chemical 
nonequilibrium [81,82], They again began with the one-dimensional Euler equations, but 
they appended species continuity equations to account for each chemical species present 
in the reacting flow and vibrational energy conservation equations to account for those 
species in vibrational nonequilibrium. The authors then redeveloped the relationships 
described previously that were required to implement Steger- Warming, van Leer, and Roe 
flux splitting. Once these splitting approaches had been implemented, a finite volume 
scheme was used along with either an explicit Runge-Kutta time integration or an implicit 
Euler time integration to solve the governing equations. The shock tube problem was again 
revisited, but nonequilibrium effects were modeled with a five species, five reaction model 
that included iV 2 ,0 2 ,iV0,iV, and O. Finally, a supersonic nozzle flow undergoing Lf 2 -air 
chemistry was modeled using a five species, two reaction model involving /f 2 , 0 2 , OH , Lf 2 0, 
and iV 2 . The calculations exhibited excellent agreement with those provided by another 
program utilizing a conventional method [81]. Extensions of the algorithm to two- and 
three-dimensions are currently underway, and once in place, this approach should also offer 
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an attractive option for modeling scramjet combustor flow fields. 

A class of more exact flux split algorithms for nonequilibrium chemistry was recently 
developed by Liu and Vinokur [83]. Steger- Warming, van Leer, and Roe flux splitting 
were again generalized for the nonequilibrium case. No simplifying assumptions were made, 
however. The most general thermal and chemical nonequilibrium flow of an arbitrary gas 
was considered. The results have not yet been incorporated into a computer program, but 
one should become available in the future. 

Additional interesting work using flux splitting has also recently been completed by 
Liou, van Leer, and Shuen [84]. This work has only recently been submitted for publication 
and so details will not be given here. The authors again employed van Leer flux vector 
splitting or Roe flux difference splitting and derived real gas forms of these approaches. 
The derivations were begun by assuming a general equation of state for a real gas in 
equilibrium. Approaches similar to those discussed previously were then used to modify 
the splittings, but the number of assumptions and particularizations employed were kept 
to a minimum. The modified splittings were then incorporated into an available TVD 
algorithm [85] and used to model several problems described by the one-dimensional 
Euler equations. Excellent agreement with the exact solution for a classical shock tube 
problem was obtained [84]. Work on the method is now continuing. 

A considerable amount of work has also been carried out over the last several years 
to develop new TVD schemes for chemically reacting real gas flows. Beginning in 1985, 
Yee developed a symmetric TVD scheme that could be employed in the context of either 
explicit or implicit numerical integration procedures [86], The approach was later gen- 
eralized to consider chemically reacting flows [87]. Yee noted that her approach could 
readily be added to existing algorithms that did not exhibit TVD behavior, e.g. the 1969 
MacCormack method, resulting in a more robust method with better shock capturing 
qualities. New explicit, semi- implicit, and implicit algorithms employing the symmetric 

TVD method were then developed and discussed [87]. An explicit multistep TVD scheme 
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was constructed using the 1969 MacCormack method [16] for the first two (predictor- 
corrector) steps followed by the addition of a conservative dissipation term as a third step, 
such that the overall scheme was TVD. The dissipative term was made up of products of 
eigenvectors of Jacobians of the governing equation system and their associated eigenval- 
ues, an entropy correction, and a limiter function. Details regarding the construction of 
the dissipative term and the determination of its magnitude were given in reference [88]. 
For the situation described earlier, where stiff kinetics resulted in small chemistry time 
scales, the explicit procedure was altered to include an implicit source term while retain- 
ing explicit flux derivatives and the same dissipative terms. Finally, a fully implicit TVD 
method was developed that included both implicit source and flux terms for situations 
where both chemistry and fluid scales were small. Calculations using both the explicit and 
semi-implicit schemes were also given in reference [87], The explicit scheme was used to 
model a shock interacting with an obstacle in a two-dimensional flow. The shock was cap- 
tured quite crisply, and more subtle flow discontinuities, including slipstreams that were 
present, were also captured. The semi-implicit scheme was generalized to three-dimensions 
and incorporated into the Uenishi semi-implicit code [56] that was described earlier in this 
chapter [89]. That program also contained the two step Rogers-Chinitz if 2 -air chemistry 
model discussed before. The three-dimensional TVD code was then used to model the 
reaction of a stoichiometric H 2 -& ir mixture flowing supersonically in a channel [89]. The 
flow was ignited by a shock produced by a wedge along the lower wall of the channel. The 
resulting fluid and chemistry profiles predicted by the program in the neighborhood of 
the shock did not exhibit the overshoots and undershoots typical of those observed with 
classical shock capturing methods. 

When stiff terms are present in the governing equation system and the fully implicit 
TVD procedure of [87] is used to solve two- or three-dimensional problems, the proce- 
dure employed to temporally integrate the governing equations must be carefully chosen. 
When ADI procedures are used, the factorization error that results when the implicit 
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operator is spatially factored often cannot be neglected [87]. Iterative procedures such 
as line Gauss-Seidel or Newton iteration, applied in the approaches discussed earlier, are 
attractive options. An alternate procedure developed by Gnoffo employed point implicit 
relaxation [90,91]. That procedure was used by Gnoffo in his three-dimensional finite 
volume code that employed a symmetric TVD upwind discretization of the governing 
Navier-Stokes, species continuity, vibrational, and electron energy equations. Pseudo-time 
relaxation was used to drive the solution to a steady state. A point-implicit procedure 
implies that during a relaxation step, the equations at each cell in the computational 
plane are discretized by using the latest available data from neighbor cells and by implic- 
itly updating only those terms which are functions of variables at the cell center. This 
procedure has proven to be very efficient on vector computers. Two options for coupling 
the governing fluid and chemistry equations, strong and weak implicit coupling, were also 
utilized. With strong implicit coupling, the complete equation set was solved as a unit, 
an approach typical of those described earlier. Weak implicit coupling involved splitting 
the fluid and chemistry equations into two groups, and applying the point-implicit method 
to each group separately during the relaxation process. The former approach was more 
physically exact, better accounting for complex wave interactions and fluid-kinetic cou- 
pling. The latter approach allowed for the relaxation strategy and time stepping to be 
tailored to the needs of the equation set [91]. Air chemistry was modeled in the program 
using an eleven species scheme that included N,0,N 2 ,0 2 , JV0,AT + ,0 + , W 2 + ,0^,ATO + , 
and e . Further details on the chemistry model and other physical modeling are given 
in reference [92], Gnoffo’s program was validated against several experiments and has 
performed well. The code was then used to successfully model the high speed external 
flow about several configurations. Because of the code structure, it should be relatively 
straightforward to add a /f 2 ~ a l r combustion model. Therefore, Gnoffo’s program, utilizing 
a symmetric TVD upwind discretization with point-relaxation appears to be an attractive 
candidate for modeling scramjet combustor flow fields as well. 
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Another attractive option to an ADI integration scheme for solving the spatially dis- 
cretized governing equations is an LU scheme that approximately splits the implicit op- 
erator into an upper and a lower operator which is independent of the dimensionality of 
the problem. Shuen and Yoon developed a scheme for solving the two-dimensional Navier- 
Stokes and species continuity equations governing chemically reacting flows that employed 
an implicit finite volume time marching LU method [93]. Details of the derivation of the 
LU scheme are given in the reference. The approach was quite attractive because, even 
though the method was fully implicit, it required only scalar diagonal inversion for solution 
of the flow equations and diagonal block inversion of the species equations. The authors 
stated that as a result, the scheme exhibited a fast convergence rate while requiring only 
about the same amount of work as that of an explicit method [93]. This advantage can be 
particularly important when problems with a large number of chemical species are being 
solved. Following development of the LU code RPLUS, an eight species, fourteen reac- 
tion chemistry model and the Baldwin-Lomax algebraic turbulence model was added to 
the program. The code was then compared with experimental data and calculations from 
other combustion programs, and it showed very good agreement with those data and com- 
putations [93]. Encouraged by their success, Yu and Shuen then extended the LU code 
to three-dimensions (RPLUS3D) [94]. A finite rate chemistry model and a turbulence 
model are currently being added to the program. Once completed, the three-dimensional 
LU code will provide another valuable tool for modeling scramjet combustor flow fields. 

The alternative methods described above for modeling reacting flows (with the ex- 
ception of Flux Corrected Transport) have generally exhibited second-order numerical 
accuracy in both space and time. Two high-order accurate methods have recently been 
developed and applied to high speed combustion problems. These methods offer improved 
accuracy as compared to lower-order methods on a given computational grid and reduced 
phase error. One method was developed by Carpenter using a fourth-order compact fi- 
nite difference scheme [59]. The scheme was initially developed by Abarbanel to accu- 
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rately solve the Euler equations in two- and three-dimensions [95]. Carpenter extended 
these ideas to the Navier-Stokes equations and used them to alter the 1969 MacCormack 
method, producing a fourth-order “compact MacCormack” scheme. The modifications 
did not change the basic structure of the MacCormack scheme, allowing it to be easily 
incorporated into existing codes using the 1969 algorithm. The modification significantly 
improved the accuracy of the algorithm, while markedly reducing the phase error. As 
a result, the improved scheme was able to crisply capture strong shocks with very little 
of the pre- and post-shock oscillations present in the old scheme. The algorithm in fact 
exhibited a TVD like behavior when capturing waves. Because of its attractiveness, the 
scheme was added to the two- and three-dimensional SPARK codes and applied to several 
supersonic chemically reacting flow problems [59]. Very good agreement with available 
experimental data was achieved, and the 3-D program was then successfully applied to a 
three-dimensional scramjet combustor. Results from that study will be described later in 
this chapter. 

High-order accurate spectral methods have also been applied to supersonic reacting 
flows. Drummond extended a Chebyshev spectral method developed to study transitioning 
flows [96,97] to include finite rate chemical reactions [98]. Spectral methods are based 
on the representation of the solution of a problem by a finite series of global functions, 
in this case Chebyshev polynomials. To apply this method to the Navier-Stokes and 
species continuity equations, the flux terms in these equations were restated in terms of 
Chebyshev series, and then the required spatial derivatives were taken. The resulting 
ordinary differential equations were then integrated with respect to time using a Runge- 
Kutta time stepping scheme. Drummond initially developed this technique for the one- 
dimensional Euler equations and species continuity equations [98]. The method was then 
extended to two-dimensions in reference [99] where a hybrid spectral-finite difference 
algorithm was used to model two-dimensional supersonic reacting flows. 

Now that both conventional and new alternate methods for modeling supersonic react- 


32 



ing flows have been examined, we will now explore the results from some of the problems 
where these methods have been applied. This should make more clear how far we have 
been able to carry these techniques to solve some of the difficult problems that we face in 
the course of designing and developing a scramjet engine. 

Applications 

Many of the computer programs described earlier in this chapter have been applied to 
study supersonic reacting flow fields. In this section, we will review some of those appli- 
cations. Space does not allow a review of each effort. Rather, a representative sample of 
the research will be provided to show the progress that has been made. We will begin 
with some basic two-dimensional studies performed in the late 1970’s. We will then ex- 
amine some of the first attempts to model simple scramjet flow fields and see how those 
studies affected later design activities. The early scramjet studies raised several basic is- 
sues regarding mechanisms controlling the mixing and chemical reaction of fuels and air 
at supersonic speeds. These issues spawned several research efforts undertaken to bet- 
ter understand high-speed combustion physics and to use that understanding to improve 
engine performance. Therefore, we will next examine the results of some of those ef- 
forts. With advancing computer resources in the mid 1980’s, scramjet calculations moved 
ahead to consider three-dimensional reacting flows. Those calculations began with simple 
three-dimensional configurations that served as model problems for the more complicated 
configurations to be later considered. We will complete our review of the applications 
by tracking the progress of those three-dimensional calculations, ending with a complex 
combustor configuration not greatly different from those being considered for advanced 
propulsion systems to be employed early in the next century. 

Early supersonic reacting flow simulations centered on configurations defined by ex- 
periments that were being performed at that time. Evans, Schexnayder, and Beach used 
the CHARNAL program as modified by Evans [14] to model a number of experiments 
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available in the late 1970’s and compared their simulations with the available data [100]. 
One case, the Burrows and Kurkov experiment [101], is shown in figure 1. Hydrogen 
was injected at Mach 1.0 through a slot in the lower wall of the configuration. Hot air 
passed by the slot at Mach 2.44. The air temperature was sufficient to ignite the hydrogen 
gas. Figure 2 shows a comparison between the computed results and the data includ- 
ing pitot pressure and reactant and product mass flow 35.6 cm downstream of the slot. 
The calculations were made using either a complete (infinitely fast) or a finite rate (eddy 
breakup) model [100], The agreement between the data and calculation is quite good for 
both pitot pressure and chemical species distributions. Encouraged by their success, the 

authors went on to examine several other experimental cases, each of which is discussed 
in reference [100]. 

While basic combustion studies continued, some of the effort was shifted to begin 
considering scramjet combustor flow fields. The first efforts considered only a part of the 
combustor, and they were still limited to two-dimensions. One example is given in figure 3 
which shows a sketch of a slot injection experiment undertaken to simulate the transverse 
injection of fuel into a scramjet combustor. In this case, helium was injected at Mach 1.0 
into a Mach 2.9 air cross-stream in a duct. Helium was used to represent gaseous hydrogen 
fuel because the test facility in which the experiment was performed could not support 
combustion. Weidner and Drummond modeled the flow field in the channel using the 
TWODLE code [22] and obtained the results shown in figure 3 [23]. Comparisons of the 
calculations with measured static pressures were fair and the comparison with helium mass 
fraction was quite good. (Interest in the slot injector case continued through the 1980’s. 
Shuen and Yoon revisited this case in 1988 using their new LU scheme in the RPLUS code 
and showed improved agreement with the static pressure data [93]. They also went on 
to examine a reacting hydrogen-air transverse jet case and reported those results in that 
reference.) Following the nonreacting study, Weidner and Drummond went on to consider 
several two-dimensional transverse injector configurations including the “staged” injector 
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configuration shown in figure 4 [23]. Staged injection of hydrogen fuel provided a means 
for producing a large region of separated flow between the injectors, thereby providing 
improved flameholding in the combustor. The final configuration resulting from the study 
is shown in figure 5. A large region of separation was produced between the injectors and 
significant reaction of hydrogen fuel and air took place in that region. 

Following completion of the fuel injector studies, scramjet calculations were extended 
in 1981 to consider a model problem for a complete engine module including both inlet 
and combustor. The calculations were still constrained to two-dimensions due to available 
computer resources at that time. Drummond and Weidner used the TWODLE program 
to simulate the flow through the module shown in figure 6 [22]. The left portion of 
the module provided inlet compression and the combustor made up the right portion of 
the module. A single fuel injection strut was positioned between the two module walls. 
Gaseous hydrogen fuel (indicated by the arrows) was injected at Mach 1.1 from the strut 
walls as well as from the engine side walls. That fuel reacted with air that entered the inlet 
at Mach 5.0. The resulting velocity, static pressure and temperature contours, and water 
contours are also shown in figure 6. The air entering the inlet was turned by shocks from 
the inlet leading edges. Shocks are indicated by a coalescence of pressure contour lines 
in the figure. These shocks struck near the strut leading edge and coalesced with shocks 
produced by the strut. The resulting shocks then had sufficient strength to separate the 
boundary layer when they reached the engine sidewalls, as indicated by a reversal of the 
velocity vectors. 

The transverse hydrogen fuel injectors located downstream of the minimum can also be 
seen in the figures along with their associated flow separations leading and trailing the in- 
jectors. The flow becomes subsonic near the fuel injectors due to air flow blockage and heat 
release from chemical reaction. Some reaction takes place in the separated regions ahead 
of the injectors. Significant reaction occurs downstream of the injectors. The temperature 
rise and water production associated with the reaction can also be seen in the figures. The 
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results of the simulation show that the injectors are not of sufficient strength to penetrate 
across the flow at the point of injection. The mixing layers do not meet until around 11 
cm downstream of the injectors. In an attempt to improve the level of fuel-air mixing, 
the injectors were moved a small distance upstream of their previous location. It was 
hoped that the change would improve the amount of injector interaction and enhance mix- 
ing. The resulting velocity field is shown in figure 7. An inlet-combustor interaction and 
choking then occured, resulting in large regions of separated flow in the inlet, and clearly, 
an unacceptable design. To examine the source of choking, the calculation was repeated 
Without chemical reaction. The calculation then proceeded normally without any signifi- 
cant separation. Therefore, it was concluded that thermal choking produced by chemical 
reaction and subsequent heat release near the injectors produced the inlet-combustor to- 
leration. While performed for a two-dimensional engine model problem, these scramjet 
simulations did demonstrate in 1981 the potential of numerical modeling for studying and 
better understanding engine flow fields. As computer resources improved, calculations 
of this type in three-dimensions were also attempted. The calculations also raised some 
basic questions regarding the very complex physics of mixing and reaction at supersonic 
speeds in a scramjet combustor. These questions resulted in several research efforts aimed 
at achieving an improved understanding of these phenomena. We will now move on to 
describe some of those basic studies and then complete this section by examining some 
three-dimensional simulations of scramjet combustor flow fields. 

The supersonic velocities that must exist in a scramjet combustor produce a perplexing, 
but interesting problem for the engine designer. At such high speeds, the mixing of the 
hydrogen fuel injected into the combustor and air from the inlet is reduced relative to the 
mixing achieved in lower speed operation. As a consequence of the reduced mixing, the 
degree of reaction that can occur and the overall combustor efficiency is also suppressed. 
The phenomena of reduced mixing in nonreacting mixing layer flows at supersonic speeds 
had been observed experimentally by Brown and Roshko as early as 1974 [102], and the 
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effect was further studied by Papamoschou and Roshko [103] in 1986. Numerical stud- 
ies of the nonreacting problem by Oh in 1974 [104] and Hussaini, Collier, and Bushnell 
in 1986 [105] were performed to better understand this phenomena. Oh suggested that 
the reduced mixing was related to weak shocks that formed when the local Mach number 
exceeded one about vortical structures that developed in an evolving mixing layer. He 
argued that other vortices would then interact with these shocks and produce fluctuations 
in the flow field that could ultimately reduce the turbulent kinetic energy and the degree 
of mixing. Hussaini et al, examined the detailed interaction of a vortex convecting subson- 
ically relative to a locally supersonic flow. They showed that a transient shock structure 
which they termed an “eddy shocklet” formed as the eddy accelerated, and the shocklet 
tended to deform the eddy. Finally, do to this interation, a vortex of opposite circulation 
formed, and the length scale of the original vortex was reduced. Based on those results, 
the authors then concluded that eddy shocklets reduced turbulent mixing through both 
the production of counter fluctuating vorticity and a reduction of turbulence scale. 

The effects of chemical reaction on mixing in a supersonic flow have also been con- 
sidered. Beginning in 1986, Drummond studied the supersonic reacting mixing layer that 
formed between coflowing streams of hydrogen gas and air at moderate Mach numbers in 
the range of Mach 2 [52,99,106]. His simulations using the SPARK code indicated that 
supersonic reacting flows exhibited some of the same features observed for subsonic react- 
ing and nonreacting flows. Vortical structure, noted in much of the subsonic nonreacting 
flow literature, was shown to be quite predominant. In agreement with the earlier react- 
ing subsonic literature, the vortical structure had a marked effect on chemical reaction in 
supersonic flow. Significant burning took palce in the eddies on the edges of the mixing 
layer, broadening the reaction zone relative to the layer thickness defined by the velocity 
gradient. In addition, the vortical flow resulted in the roll up of unburned reactants inside a 
layer of partially or fully burned products. This phenomena, often termed “unmixedness” 
in subsonic flows, prohibited the reaction of captured reactants and reduced the overall 
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efficiency of the combustion process. There was also some indication that heat release re- 
sulting from the chemical reaction also reduced the amount of mixing and reaction further 
downstream in a mixing layer. 

Work is continuing to better understand the phenomena that produce reduced com- 
bustion efficiency in supersonic combustors. Related efforts also began in 1987 to develop 
techniques for enhancing the degree of fuel-air mixing and combustion that occured in high 
Mach number flows. Guirguis used the flux-corrected transport algorithm of Boris to study 
convective mixing in high-speed nonreacting flows [108]. He again chose the mixing layer 
as a model problem for mixing in a scramjet combustor and solved the two-dimensional 
Euler equations to examine the effects of imposed pressure gradients on the spatial de- 
velopment of the layer. The upper stream of the mixing layer, made up of “species 1” 
entered a confined channel at Mach 4.5 and a pressure of 4 atm. The lower stream of 
species 2 entered at Mach 1.5 and a pressure of 2 atm. The stream temperatures were 
matched. In addition, two other cases were considered. In the second case, the pressures 
of the two streams were matched and the mixing layer was unconfined, and in the third 
case, the lower stream entered at a pressure of 2.1 times that of the upper stream. Other 
minor differences between the cases were discussed in reference [108]. One result from the 
study by Guirguis is shown in figure 8. In that figure, mass fraction contours of species 1 
from 0.01 to 0.99 are plotted for the three cases, respectively. The matched pressure case 
showed little development with increasing streamwise coordinate, whereas the first and 
third cases with the imposed transverse gradient showed significant spread of the layer. 
Guirguis therefore concluded that in a channel of a given length, differences in pressure 
across the layer enhanced mixing. He also suggested means whereby pressure differences 
could be induced by the geometry of the combustor. Guirguis then continued his work to 
find a means of enhancing mixing where the two stream pressures were matched [109]. He 
concluded that mixing could be enhanced if a bluff body was placed at the trailing edge of 
the splitter plate separating the two matched pressure streams. One result from the bluff 
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body case is shown in figure 9. The mass fraction contours given in the figure again show 
a significant degree of enhancement is provided by the presence of the bluff body. 

Kumar, Bushnell, and Hussaini also examined the problem of mixing in flows undergo- 
ing supersonic combustion using a two dimensional version of the NASCRIN Code [110]. 
Several techniques for enhancing turbulence and mixing were suggested and one enhance- 
ment technique that employed an oscillating shock was studied numerically. In this case a 
premixed stoichiometric hydrogen-air flow was processed through a spatially and tempo- 
rally oscillating shock wave, and the resulting flow was studied with and without chemical 
reaction. Reaction was modeled using an equilibrium chemistry model. One reacting flow 
case that was considered is shown schematically in figure 10. Here, the premixed fuel-air 
flow was introduced into a two-dimensional channel at Mach 3, a pressure of 1 atm, and a 
temperature of 1500 K, and the mixture was passed through a 10 degree shock produced 
by a lower wall compression. A periodic oscillation was imposed upon the flow near the 
lower wall, and this resulted in a periodic oscillation of the shock wave. The resulting 
pressure field over one period of the imposed oscillation is shown in figure 11. A wave can 
be seen to propagate along the shock. The oscillating shock was shown to increase the 
level of turbulence in the flow field, and the degree of turbulence enhancement was seen to 
increase with a decreasing frequency of shock oscillation. Chemical reaction as defined by 
the equilibrium model was shown to have little effect relative to the nonreacting results. 

Several techniques for enhancing mixing and reaction in a model scramjet combustor 
were studied in 1988 by Drummond and Mukunda [107]. They considered the shock and 
expansion structure that existed in scramjet combustors, and sought to redirect this exist- 
ing wave structure to utilize it for enhancement. Planar and curved shocks were created by 
structure and fuel injection in scramjet engines, and so both shock shapes were considered 
as candidates for mixing enhancement. Three cases that were studied to determine the 
efficiency of shock shapes are shown in figure 12. Case 1 shows a two-dimensional spatially 
developing supersonic mixing layer that served as a baseline to indicate the degree of mix- 


39 



ing and reaction that could be achieved without enhancement. Here, gaseous hydrogen 
fuel and air are initially separated by a thin splitter plate. Hydrogen and air both enter 
at Mach 2, a pressure of 1 atm, and a temperature of 2000 K. The hydrogen and air begin 
mixing at the trailing edge of the splitter plate, ignition occurs at some small distance 
downstream, and combustion continues from that point. Two enhancement studies, cases 
2 and 3, are also shown in figure 12. In case 2, the mixing layer is processed through 
two 10 degree planar shocks produced by wedges in the fuel and air streams. In case 3, 
a small interference body that produces a curved bow shock is placed at the center of 
the mixing layer. The flow field in each case was then simulated using the SPARK code. 
Chemistry was modeled using either a 4 species, 1 reaction or a 9 species, 18 reaction 
H 2 -a\r finite-rate chemistry model. The resulting instantaneous mixing efficiencies shown 
as a function of streamwise coordinate are given in figure 13 [107]. The degree of mixing 
achieved in the mixing layer without enhancement is shown by the result from case 1. 
Very little additional mixing results from passing the layer through planar shocks in case 
2. Case 3 exhibits a significantly higher degree of mixing. When the high velocity gradient 
across the mixing layer is processed by the curved bow shock, vorticity is produced, and 
the mixing is subsequently enhanced downstream of the shock. 

With a background of planar and curved shock mixing enhancement in hand, fuel 
injection configurations typical of those used in a scramjet engine were considered. A 
typical design is described in figure 14, which shows the trailing edge of a NASA Langley 
fuel injection strut. Inlet air crossed the strut as shown. Gaseous hydrogen fuel was 
injected parallel to inlet air from the base of the strut, and transverse to the inlet air 
behind a small rearward facing step in the surface of the strut. (The step caused a small 
recirculation region to form which captured a small amount of fuel and air. Combustion 
then occured in that region providing a flameholding sight in the engine combustor.) The 
results of the previous analysis suggested a simple change in the fuel injector configuration 
that might provide improved fuel-air mixing and enhanced combustion efficiency. If the 
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parallel injector from the strut base were moved to the step face as shown in figure 14, 
then the high velocity gradient of that jet would interact with the curved bow shock that 
formed ahead of the transverse injector. The mixing enhancement modification was then 
examined by again simulating the two-dimensional flow fields. Chemical reaction was 
described by the 9 species, 18 reaction finite rate model. The resulting water mass fraction 
contours downstream of the modified strut for only the parallel injector and for both the 
parallel and transverse injector are given in figure 14. The multiple jet configuration clearly 
gave significantly higher levels of jet development, mixing, and reaction. A quantitative 
assessment of mixing and reaction is given in figure 15. In that figure, the combustion 
efficiency for both cases is plotted against the streamwise coordinate. Jet interaction 
produced a significantly more rapid increase in mixing and combustion. Instantaneous 
combustion efficiencies of about 80 percent were obtained for the enhanced case within 
the first 40 percent of the solution length whereas the unenhanced case required about 
75 percent of the solution length to achieve the same level of reaction. Enhancement 
by such approaches could therefore result shorter combustor lengths for a given level of 
performance. 

While the basic studies we have discussed were underway, several interesting three- 
dimensional scramjet combustor simulations were undertaken in 1987 and 1988. Uenishi, 
Rogers, and Northam used their three-dimensional combustor code, described earlier, to 
study several generic combustor configurations that utilized wall fuel injection [56,58]. One 
of the interesting cases is shown in figure 16. This combustor model was also being used in 
an experimental program at NASA Langley when the study was undertaken, although no 
data was yet available. Vitiated hot air, produced by the combustion of hydrogen, oxygen, 
and air in a facility heater, entered the model combustor from the left. Gaseous hydrogen 
fuel was initially injected a small distance downstream from this point through a row of 
secondary orifice fuel injectors located on opposite walls of the combustor. These injectors 
were intended to pilot a pair of primary hydrogen fuel injectors again located on opposite 
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walls some distance downstream of rearward facing steps in the walls. The flow in the 
combustor was modeled using the Uenishi code [58]. Features of the combustor flow field 
resulting from that study are shown in figure 17. Figure 17a shows the distribution of 
the injected fuel given as the mass fraction of total H 2 in the streamwise x-z plane of the 
computation along the duct centerline. The total H 2 distribution represented the sum at 
each point of the hydrogen atoms in any form, i.e. H 2} H 2 0 , or OH , so that it represented 
the amount of injected H 2 fuel. The fuel rich regions in the vicinity of the injectors are 
clearly evident as is the penetration of the fuel across the duct. The contours in figure 
17a show that the upstream fuel injectors provide relatively well mixed fuel and air into 
the downstream region which then arguments the fuel mixing and combustion process 
produced by the downstream injectors [58]. 

The three-dimensionality of the flow field is clearly shown by the distributions of total 
H 2 and static temperature given in figure 17b at several cross-planes along the combustor 
length. The static temperature was normalized by the initial static temperature. The 
locations of the y-z cross-planes are indicated in figure 17a. The distribution of total H 2 
in each cross-plane indicates the mixing of injected fuel. At location A, the fuel jets from 
the upstream injectors have not merged, and the cold fuel has not mixed with the hot air 
stream to an extent sufficient to result in any significant heat release, as is indicated by the 
absence of any appreciable temperature rise at this location. At location B, however, the 
fuel has mixed and reacted with the hot-air stream to an extent that the temperature is 
about 1.7 times the air stream. Further downstream, at locations C and D, the 0.5 percent 
contour line becomes more uniform and a large area of the flow contains reacting fuel as 
also indicated by the elevated temperatures. Results from location E indicate that some 
fuel from the downstream injector is entrained upstream by the recirculating flow behind 
the step. This upstream entrainment of H 2 from the downstream injector is also indicated 
by the regions of high and low temperature at location E. The fuel rich regions have a 
lower temperature because the H 2 is cold and because the fuel-air mixture is too rich to 
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react. At locations F and G, which are downstream of the primary fuel injector, large 
regions of unreacted fuel can be seen. With the present flow conditions, the H 2 injected 

from the primary injectors penetrates very well, passing beyond the combustor centerline 
[58]. 

Calculations are continuing with the Uenishi code to model complex three-dimensional 
combustor flow fields. A major effort is now also underway to compare the code with addi- 
tional experimental data that has recently become available. Work is also currently being 
completed to couple the spatially elliptic Uenishi code to several parabolized Navier-Stokes 
programs, including the Chitsomboon PNS code and the Kamath PNS code, that were 
described earlier in this chapter. Parabolized Navier-Stokes codes can be more efficient for 
simulating the combustor far-field well away from the fuel injectors where the flow is spa- 
tially elliptic due to separation. Therefore, there is an incentive, based on computational 
efficiency, to switch to PNS programs wherever possible in the combustor. 

A number of interesting combustor calculations have been made by Gielda using his 
three-dimensional parabolized Navier-Stokes program, SSCPNS, that was also described 
earlier in this chapter. One configuration that he analyzed is given in figure 18, which 
shows a sketch of a generic hypersonic propulsion system including an inlet, combustor, 
and nozzle [67]. Conditions entering the inlet were determined by assuming that the 
hypersonic vehicle was cruising at Mach 16 and at an altitude resulting in a free stream 
temperature and pressure of 452 R and 0.004 atm, respectively. Gaseous H 2 fuel was 
injected through five injectors into the engine beginning at a location x/L = 0.625 and 
continuing for a distance of x/L = 0.24 along the combustor as indicated by figure 18. 
This injection schedule resulted in a fuel equivalence ratio of approximately three. The 
fuel was injected at an angle of 25 degrees to the engine cowl defined in the figure by the 
x-axis. 

Once the engine geometry and fuel injection configuration had been defined, the entire 

engine flow field was simulated using the SSCPNS code. The computed flow field was 
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highly three-dimensional. The three-dimensionality can be seen in figure 19 which shows 
Mach number contours at a number of cross-planes along the length of the inlet, combustor, 
and enclosed portion of the nozzle. Gielda’s study also showed that there were significant 
losses associated with the fuel injection process in the generic engine. The losses resulted 
from the formation of strong oblique shocks as the inlet air interacted with the injected 
fuel. The induced shock waves can be seen in figure 20 which shows the computed pressure 
contours m the engine. The entropy increase that resulted when the engine flow passed 
through these oblique shocks caused a significant loss in the average stagnation pressure 
of the flow. This loss in stagnation pressure would ultimately result in a loss in the overall 
thrust that could be achieved. The fuel injector configuration did result, however, in good 
fuel-air mixing and combustion. The degree of combustion is indicated by the computed 
water mass fraction contours shown in figure 21 at several cross-planes along the combustor 
and nozzle lengths. Water mass fractions of up to 0.23 resulted across a large extent of the 
cross-planes, and an oxygen utilization of approximately 0.80 was achieved. Normalized 
static temperature contours in the engine resulting from the simulation are also shown in 
figure 21 [67]. 

In order to improve the performance of the engine, Gielda further investigated the 
fuel injection process in an attempt lower the associated shock losses. Three numerical 
experiments were performed that included a study of the engine flow field without injection 
or reaction, with injection, and with injection and reaction. By comparing the three cases, 
Gielda was able to assess the losses associated with the fuel injection process and associated 
with the heat addition resulting from the combustion process [67]. He concluded that 
further investigation into the injection process was required and that such a study would 
require simulation of the flow field near the fuel injectors using the spatially elliptic form 
of the Navier-Stokes equations. Such a study has recently been undertaken. 

Carpenter used the SPARK3D spatially elliptic Navier-Stokes program, described ear- 
lier, to model the generic scramjet combustor sketched in figure 22 [59]. This combustor 
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configuration has been used as a generic model for actual combustor designs being consid- 
ered as a part of a propulsion system for a hypersonic vehicle that will become operational 
early in the next century. Flow processed by an inlet enters the combustor from the left in 
the figure. The flow exits the inlet at a velocity of 1500 m/s, a static temperature of 1000 
K, and a pressure of 0.5 atm. That flow then passes over rearward facing steps on the walls, 
experiencing a mild area expansion. Gaseous H 2 fuel is injected from a line of four orifice 
fuel injectors that lie along the face of the step on each of the four engine walls. The fuel 
and inlet air begin mixing just downstream of the steps. The flow for this configuration 
was simulated by Carpenter using the SPARK3D code [59]. Because the combustor was 
symmetric about its center, only the lower quadrant, shown in figure 22 was modeled. In 
addition, even though the program contained a generalized finite rate chemistry model, 
only fuel-air mixing was considered at this stage of the study. The resulting H 2 mass 
fraction distributions at 2 cm and 5 cm downstream of the step are given in figure 23. At 
the first station, recirculation behind the step and diffusion spread the hydrogen towards 
the walls, but there is little penetration of the fuel into the core flow in the center of the 
combustor. The 0.2 mass fraction contour does not extend beyond the height of the step. 
(Fuel rich regions having contours with a mass fraction greater than 0.2 were surpressed 
for clarity.) Further downstream at the 5 cm station, hydrogen mixes to a greater extent 
with the air, and the mixture extends further into the core flow. The extent of mixing 
is still not sufficient to provide an acceptable level of mixing efficiency or for that matter 
combustion efficiency if reaction had been allowed. As a result of these findings, mixing 
enhancement strategies similar to those discussed earlier are currently being applied to the 
generic combustor configuration to achieve an acceptable level of combustion efficiency. 

Several other interesting calculations have also recently been undertaken to model 
scramjet combustor flow fields. But we have now reached the end of 1988, and it is time to 
bring this chapter to a conclusion. Computational methods have clearly been established 
as valuable tools for simulating combustor flows and providing an improved understanding 
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of their complex physics. Therefore, this story will continue to be told! With expected 
improvements in algorithms, computer resources, and creativity, the story should only 
improve, and computational methods should take on an even greater role in studying and 
designing high speed propulsion systems. 

Concluding Remarks 

Significant progress has been made in the past twenty years to develop techniques for 
modeling supersonic reacting flows. The progress has been due both to an advancement 
in numerical methods for computing reacting flow fields as well as an appreciable growth 
in computer speed and storage. The improvements in both algorithms and computer 
power have also lead to a gradually improved understanding of the physics of reacting 
flows which 1 8 so very important if computation is to have a truly significant impact on 
scramjet combustor design. The national program to develop a trans-atmospheric vehicle 
has been quite effective in generating a renewed interest in supersonic combustion, and 
it has spawned several programs aimed at better understanding its complexities. While 
these programs have begun to unravel some of the mysteries, they have also shown us that 
supersonic reacting flows still contain many unsolved questions. These questions provide 
exciting opportunities for future research, and these opportunities should not be missed! 

The efficient mixing and combustion of fuel and air in a scramjet combustor remains 
a critical issue today. We are just beginning to understand the phenomena that suppress 
mixing in these high Mach number devices, but a great deal of work remains to be done. 
Once the mixing question is answered, many opportunities will still remain to employ this 
improved physical understanding to produce an efficient mixing strategy. Ground based 
experimental facilities will be used to study and design scramjet combustors by providing 
flows consistent with those that the engine will injest over its operating envelope. Ground 
based facilities are only able to create continuous flow conditions up to a flight Mach num- 
ber of about 8, however. Beyond Mach 8, the only options available to the experimentalist 
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are pulse facilities that create flight conditions for only a short period of time, or the very 
vehicle that the national program seeks to build! Numerical methods provide an alter- 
native to the Mach 8 barrier, but only if they are properly applied to the problem. To 
apply these methods properly, a fundamental understanding of the physical phenomena 
present in high speed reacting flows must be achieved. Even with this understanding, we 
will still be required to simulate flow fields in large geometric configurations, necessitating 
highly efficient numerical algorithms and significant advancements in computer technology. 
There is certainly a problem here for everyone! The problems are difficult ones, but they 
can be solved if we continue to tackle them in a careful and dedicated way. The United 
States committed itself to landing men on the moon in 1961, and it succeeded through a 
dedicated national commitment to do so in only eight short years. The national program 
to build and fly a trans-atmospheric vehicle is no more difficult a goal to achieve. It simply 
requires another national commitment to succeed, and a willingness on our part to accept 
the exciting challenges that such a program provides. 
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TABLE 1. - Finite-Rate Chemistry Model and Arrhenius Rate Coefficients for Each Reaction 


Reaction 

number 

Reaction 

A 

N 

E, kJ/g-mole 

1 . 

h 2 +o 2 =oh+oh 

0.1700E+14 

0.00 

201.5 

2. 

h+o 2 =oh+o 

0.1420E+15 

0.00 

68.6 

3. 

oh+h 2 =h 2 o+h 

0.3160E+08 

1.80 

12.78 

4. 

o+h 2 =oh+h 

0.2070E+15 

0.00 

57.5 

5. 

oh+oh=h 2 o+o 

0.5500E+14 

0.00 

29.3 

6. 

h+oh=h 2 o+m 

0.2210E+23 

-2.00 

0.0 

7. 

h+h=h 2 +m 

0.6530E+18 

-1.00 

0.0 

8. 

h+o 2 =ho 2 +m 

0.3200E+19 

-1.00 

0.0 

9. 

ho 2 +oh=h 2 o+o 2 

0.5000E+14 

0.00 

4.2 

10. 

ho 2 +h=h 2 +o 2 

0.2530E+14 

0.00 

2.9 

11. 

ho 2 +h=oh+oh 

0.1990E+15 

0.00 

7.5 

12. 

ho 2 +o=oh+o 2 

0.5000E+14 

0.00 

4.2 

13. 

ho 2 +ho 2 =h 2 o 2 +o 2 

0.1990E+13 

0.00 

0.0 

14. 

ho 2 +h 2 = h 2 o 2 +h 

0.3010E+12 

0.00 

78.2 

15. 

h 2 o 2 +oh=ho 2 +h 2 o 

0.1020E+14 

0.00 

7.9 

16. 

h 2 o 2 +h=oh+h 2 o 

0.5000E+15 

0.00 

41.9 

17. 

h 2 o 2 +o=oh+ho 2 

0.1990E+14 

0.00 

24.7 

18. 

m+h 2 o 2 =oh+oh 

0.1210E+18 

0.00 

190.4 

For the single step reaction 

0.5510E+15 

0.00 

30.2 



A i 

b< 

C, 

C, 

c p 

Da 

D t 

E 

E f F t G 

h 

9j 

Gr 

h< 

h° 

H 

K b 

K f 

K eq 

M 

Mi 

n { 

ns 

nr 

P 
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Appendix 

Nomenclature 

reaction rate constant for jth reaction 

body force of species i 

concentration of species i 

time rate of change of C{ 

specific heat at constant pressure 

binary diffusion coefficient 

thermal diffusion coefficient 

total internal energy; activation energy 

flux vectors in x, y, and z coordinate directions 

mass fraction of species i 

Gibbs energy of species i 

Gibbs energy of reaction 

enthalpy of species i 

reference enthalpy of species i 

source vector 

backward rate constant 

forward rate constant 

equilibrium constant 

Mach number 

molecular weight of species i 
moles of species i 
number of chemical species 
number of chemical reactions 
pressure 
heat flux 
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R° 

T 

T r 

t € 

t 

At 

— ♦ 
U 

u 

v 

Ui 


Vi 

Vi 

Wi 

X 

y 

Xi 

X 

in 

7 

<5 

P 




T 


Pi 

P 

t 


universal gas constant 
temperature 

reference temperature = 298 °k 

effective temperature 

time 

time step 

dependent variable vector 
streamwise velocity 
transverse velocity 

streamwise diffusion velocity of species i 

transverse diffusion velocity of species i 

diffusion velocity vector of species i 

species production rate of species i 

streamwise coordinate 

transverse coordinate 

mole fraction of species i 

second viscosity coefficient 

stoichiometric coefficient; species i, reaction j 

ratio of specific heats 

Kronecker delta function 

density 

normal stress 

effective collision diameter 
shear stress 

laminar viscosity of species i 
mixture laminar viscosity 
computational streamwise coordinate 
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computational transverse coordinate 
diffusion collision integral 
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Alratraami M__ ■ 1*75 P„ - 1.21 at* T** ■ 1212 1C 

FS PS rs 

1 . - 1.5h* ft - .25 

B.L. 


(Jpatraam Injactorai 


Downatraam Injactorai 


T w.H * *50 X 


- 1.0 Pj - 2.1 atm Tj - 250 K 

A • .1 p - .57 C - .9 
r d 


- 1,0 Pj - 12.0 atm Tj - 250 K 

A . .78 p - 3.2 C. • .93 
r d 
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L = 20.0 cm, W » 10.0 cm 
S » 2.0 cm, Lj- 1.5 cm, D » 3.5 mm 

T * 1000 K, P « 0.5 atm, V • 1500 m/s 
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